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Abstract 

We present a finite difference method to solve a new type of nonlocal hydro- 
dynamic equations that arise in the theory of spatially inhomogeneous Bloch 
oscillations in semiconductor superlattices. The hydrodynamic equations de- 
scribe the evolution of the electron density, electric field and the complex 
amplitude of the Bloch oscillations for the electron current density and the 
mean energy density. These equations contain averages over the Bloch phase 
which are integrals of the unknown electric field and are derived by singular 
perturbation methods. Among the solutions of the hydrodynamic equations, 
at a 70 K lattice temperature, there are spatially inhomogeneous Bloch os- 
cillations coexisting with moving electric field domains and Gunn-type oscil- 
lations of the current. At higher temperature (300 K) only Bloch oscillations 
remain. These novel solutions are found for restitution coefficients in a nar- 
row interval below their critical values and disappear for larger values. We 
use an efficient numerical method based on an implicit second-order finite 
difference scheme for both the electric field equation (of drift- diffusion type) 
and the parabolic equation for the complex amplitude. Double integrals ap- 
pearing in the nonlocal hydrodynamic equations are calculated by means of 
expansions in modified Bessel functions. We use numerical simulations to as- 
certain the convergence of the method. If the complex amplitude equation is 



* Corresponding author. 
Email addresses: mariaiio.alvaroQuc3m.es (M. Alvaro), 
manuel.carretero@uc3m.es (M. Carretero), bonilla@ing.uc3m.es (L. L. Bonilla) 
URL: http://scala.uc3m.es (L. L. Bonilla) 



Preprint submitted to Elsevier 



March 3, 2013 



solved using a first order scheme for restitution coefficients near their critical 
values, a spurious convection arises that annihilates the complex amplitude 
in the part of the superlattice that is closer to the cathode. This numerical 
artifact disappears if the space step is appropriately reduced or we use the 
second-order numerical scheme. 

Keywords: Semiconductor superlattice, Bloch oscillations, nonlocal 
hydrodynamic equations, spurious convection, self-sustained current 
oscillations 
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1. Introduction 

An immediate consequence of the Bloch theorem is that the position of 
an electron inside an energy band of a crystal oscillates coherently under 
an applied constant electric field with a frequency proportional to the field 
To observe these so-called Bloch oscillations (BOs), their period must be 
much shorter than the scattering time for, otherwise, scattering eliminates 
them. The required electric field to observe Bloch oscillations is too large 
for a natural crystal, but it becomes reasonable if artificial crystals with 
larger spatial periods are created. An artificial crystal can be formed by 
growing many identical periods comprising a number of layers of two different 
materials. The resulting superlattice (SL) was suggested by Esaki and Tsu 
in 1970 as a possible realization of Bloch oscillations Damped Bloch 
oscillations were first observed in 1992 in semiconductor SLs whose initial 
state was prepared optically j^. In recent years, BOs have been observed 
in other artificial crystals such as atoms placed in the potential minima of 
a laser- induced optical standing wave 0, Isf, photons in a periodic array of 
waveguides 0, 0] and Bose-Einstein condensates in optical lattices [8] among 
other systems |9|. 

In SLs made out of doped semiconductors, scattering usually destroys 
BOs, but we have shown recently that BOs can persist even in the hydrody- 



namic regime for a SL with long scattering times [10|, O]. To do so, we con- 
sider a Boltzmann-Poisson description of a SL with a single occupied electron 
miniband and a dissipative Bhatnagar-Gross-Krook (BGK) collision model 



12[. In semiconductor SLs, collisions are inelastic: they conserve charge but 
dissipate energy and momentum. Then the BGK inelastic-collision model 
contains two restitution coefficients that regulate the fractions of energy and 
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momentum lost in collisions 12|. Previously, an inelastic BGK model was 



used to describe gran ular fluids, in which energy (but not momentum) is lost 



during collisions [13|. Boltzmann-BGK mass- and energy-conserving kinetic 
equations with different local equilibrium distributions describe a particle in 
an external potential, and have been used to derive energy-transport equa- 
tions and to prove H theorems; see IJ] and references cited therein. Quantum 
hydrodynamic and quantum energy-transport models have been derived from 
quantum BGK kinetic equations that conserve mass, momentum and energy; 
see [lil and references cited therein. For SLs, inelastic collisions imply that 
BOs are damped on a time scale given by the scattering time. Damping mod- 
ulates the amplitude of the BOs if the scattering time is sufficiently long. In 
the limit of almost elastic scattering and high electric field, the electron cur- 
rent density and mean energy oscillate at the Bloch frequency, whereas the 
electron density, the electric field and the envelope of the BOs vary on a 
slower scale. In this limit, it is possible to derive a set of one- dimensional 
nonlocal hydrodynamic equations for the electric field and the complex am- 
plitude of the BOs. Their solutions allow to reconstruct the rapidly varying 
electron current and mean energy densities. The hydrodynamic equations 
are of a quite novel type: they contain averages over the phase of the BOs 
which is proportional to the integral over time of the electric field, and there- 
fore the phase is also an unknown to be determined. Appropriate boundary 
and initial conditions include initiation of the BOs possibly by optical means 
[ij. Numerical solution of the hydrodynamic equations shows stable BOs for 



appropriate parameter values [10|, lUj. For dc voltage biased SLs at suffi 



ciently low lattice temperatures, there are solutions in which the amplitude 
of the BOs is time-periodic and the electric field profile inside the SL exhibits 
electric field domains (EFDs) 

In this paper, we present a method to solve numerically the nonlocal hy- 
drodynamic equations describing BOs for a dc voltage biased SL. Although 
the problem is one- dimensional (ID), it is very computationally intensive, so 
we need an efficient numerical method to solve it. We solve the hydrody- 
namics equations by means of an efficient implicit finite difference numeri- 
cal scheme, which uses a fixed point iteration process to obtain numerically 
both the electric field and the BO complex amplitude at each time step. 
The equation for the field is a nonlocal drift-diffusion equation (containing 
integrals over the Bloch phase which is an integral of the electric field) that 
is solved using an implicit numerical scheme involving the inversion of only 
one tridiagonal matrix per iteration. This equation is coupled to that of the 
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BO amplitude. We use a second order numerical scheme to solve the latter. 
In the hydrodynamic equations, there appear several Fourier coefficients of 
the Boltzmann distribution function (which is periodic in the momentum 
variable with period 27r/Hf Z is the spatial period of the SL). These Fourier 
coefficients are approximated by truncated series of modified Bessel functions 
for computational efficiency. With appropriate parameter values and bound- 
ary conditions, numerical solutions show that initial profiles for the field and 
the BO amplitude evolve to stable spatially inhomogeneous proffies at room 



temperature At low temperature (70 K), we have found that BOs and 



Gunn-type oscillations [16|,ll7[ due to EFD dynamics may coexist. Increasing 
lattice temperature produces large diffusion coefficients as compared to the 
convective part of the averaged drift-diffusion equation for the electric field. 
This eliminates the Gunn-type oscillations. At low lattice temperature, the 
diffusion does not change that much, but convection dominates the average 
electron current density, thereby facilitating movable EFDs and Gunn-type 
oscillations BOs (accompanied or not by Gunn-type oscillations in their 
amplitude) disappear as the collisions in the kinetic equation become more 
inelastic and the BO amplitude becomes zero everywhere. If the amplitude 
of the BOs is set to zero, the drift-diffusion equation for the electric field is 
similar to those obtained with a local equilibrium distribution that depends 



only on the electron density [18|. Solutions of this drift-diffusion equation 
include Gunn-type self-oscillations due to EFD dynamics. Direct solution of 
the Boltzmann- Poisson system studied in |18j con&ms this [lij. 

The rest of the paper is as follows. In Section |2l we describe the Boltzmann- 
BGK-Poisson system and the nondimensional hydrodynamic equations de- 
rived from it. In Section [3l we explain the numerical method for solving 
the hydrodynamic equations as well as the numerical results. The analysis 
of the convergence of the numerical method is based on numerical simula- 
tions and it is presented in section |H Section |5] contains our conclusions. In 



Appendix A we include some technical details including the series of modi- 



fied Bessel functions used to approximate some integrals. 

2. Model equations 

For a semiconductor superlattice with a single occupied miniband of dis- 
persion relation 

g(k) = ^{1- cos kl) (1) 
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(A is the SL miniband width and / is the SL spatial period), the distribution 
function f{x,k,t) of electrons with position in the interval (x, x + dx) and 
wave vector in {k, k + dk) satisfies the system of equations [ij, 
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Here f^, n, N^, £, v{k) = h^^d£/dk, ks, — e < 0, Ua, Ta, Jn, E, Eq 
and —F = —dW/dx are the local equilibrium distribution, the 2D electron 
density, the 2D doping density, the dielectric constant, the group velocity, 
the Boltzmann constant, the electron charge, the hydrodynamic velocity, the 
electron temperature, the electron current density, the mean energy density, 
the lattice mean energy density, and the electric field, respectively. W is the 
electric potential. The lattice mean energy density is related to the lattice 
temperature Tq as Eq = ^ /] ( 2k^To ) /-^o ( 2k^To ^ ' "^^ere Ij{x) is the modified 
Bessel function of index j [11]. Note that the ID distribution functions / 
and have the same units as the 2D electron density n and that Ua and 
are functions of n, Jn and E obtained by solving ©-([H]) with given by 
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(jni). The ID Boltzmann local equilibrium is an approximation to a more 



general Fermi-Dirac distribution u is the collision frequency which we 
take as a constant for the sake of simplicity, ae and aj are constant restitution 
coefficients that indicate the fraction of energy and momentum dissipated in 
inelastic collisions (with phonons, for example). The distribution function is 
periodic in k with period 27r//. 

Ktitorov, Simin and Sindalovskii (KSS) ji^] proposed in 1972 an equa- 
tion similar to ([2]) except that was replaced by the Boltzmann equi- 
librium distribution at the lattice temperature Tq and an additional term 
Qp = —h'p[f{x,k,t) — f{x. — k,t)]/2 (representing ID impurity collisions) was 



added to the RHS. Later Ignatov and Shashkin [2l| proposed a local equi 



librium distribution function similar to with Ua = and = Tq. Such 



a kinetic model was numerically solved by Cebrian et al [19| using a particle 
method, and Bonilla et al p^] derived from it a generalized drift- diffusion 
equation for the electric field exhibiting Gunn-like oscillations of the current 
due to EFD dynamics. However the KSS model with the Ignatov-Shashkin 
modification cannot sustain stable Bloch oscillations [lH. To see the relation 
with BOs, we multiply ([2]) by 1, v{k) or [A/2 — S{k)] and integrate over k, 
thereby obtaining the following moment equations for n, J„ and E: 

dJ eAH d eHuEF 

^ + w 9^^" - — w- = -""^-^^ ^''^ 

dE_lEdJn_^ + ^ = -vaJE - Eo) (13) 

dt en dx 8hn dx n ^ 

Here we have used ([1]) and the Fourier coefficients fj of the periodic distri- 
bution function: 

oo 

f{x,k,t)= J2 fA^,t)e'"''. (14) 

j=-oo 

Note that J„ = -eA Im/i/(2n) and E = A Re/i/(2n). ([n]) is the charge con- 
tinuity equation. For elastic collisions, ae = aj = and space-independent 
n, Jn and E, we obtain dn/dt = (thus n is constant), and f[T^ and f[T^ 
become 

dJn eHnEF 

^-^^ = °' 
dE FJJ 

— + = 0. 16 

at n 
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Since n is constant, J„ and E are time periodic and oscillate with the Bloch 
frequency ub = eFl/h, proportional to the electric field. In the general case 
of space dependent moments, f|TTl) - f|T3|) are not a closed system of equations 
because they depend on the second harmonic of the distribution function /2. 
In an appropriate limit, the terms on the RHS of f|T2|) and f|T3|) modulate 
the BOs, so that n, F and the amplitude of the BOs evolve on a slower 



time scale. Based on these ideas, we have derived in a system of slowly- 
varying nonlocal hydrodynamic equations for these magnitudes using singular 
perturbation methods. 
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Table 1: Hyperbolic scaling and nondimensionalization with v = Hz. 



2.1. Hydrodynamic equations 

Written in the nondimensional units of table [H the hydrodynamic equa- 
tions are 11 
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F{x,s)ds. (22) 



Here 7ej- = C(e,j/S are rescaled restitution coefficients and ff and are the 
Fourier coefficients, 

ff = J' f e-"^ dk, (23) 

of the Boltzmann local equilibrium distribution ([5]) which, in nondimensional 
form, is: 

^ ^uk+P cos fc 

The Boltzmann local equilibrium distribution can be expanded in a power 
series of the small parameter 6: 

oo 
m=0 

In fl2^ . M and /3 are nondimensional multipliers corresponding to Ua and to 
l/{kBTa) in (E]). They are found by solving the nondimensional versions of 
© and dH]): 

^ [ f sink dk = {1-6^,) Jn, (27) 
f''coskdk = E-6^,{E-Eo), (28) 



when is given by fl25l) . In f|T71) and fl2Tl) . /j^^^"* is given by 

/fW=Ae-*^ (29) 
To calculate f2-i in (121]), we need m*^^^ and /S^-^^ given by 

/f « = 7eni?o + /i,5 - A e-^' - Ae^^ (30) 

where A is the complex conjugate of A, and 

B(0) 



(31) 



Together with the Poisson equation (!20|) . (ITTI) is a drift-diffusion equation 
for the electric field —F, ( !T8|) gives the time evolution of the BO complex 
amplitude A, and (fT9|) gives the total current density J (proportional to the 
electric current in the circuit attached to the superlattice) . The restitution 
coefficients are rescaled as = Sje, «j = Sjj, where S = l/(z/[t]) ^ 1 is the 
ratio between the Bloch period and the convective time scal e \t ] in table [1] 
The hydrodynamic equations hold in the limit as 5 — )■ 0+ There are 

two time scales in the hydrodynamic equations: a nonlinear fast time scale 
9 (the phase of the BOs, on the picosecond time scale), given by f l2^ . and 
a slow time scale t on the nanosecond scale. The electric field —F{x,t), the 
electron density n{x, t), the BOs complex amplitude envelope A{x, t) and the 
^-averaged current density (J)g (t) all vary on the slow time scale and may 
exhibit Gunn-type oscillations with frequencies on the GHz scale. Once we 
have obtained F, {J)g and A from (ITTI) and (fTS]) . which depend only on the 
slow time scale t, we can calculate from f ll9p the total current density J, 
which depends on both time scales. Although A = is an exact solution of 
f lTSj) . we are interested in finding numerical solutions with undamped BOs 
(i.e. A ^ 0) coexisting with Gunn type oscillations of the current. The 
first term in the RHS of f|T8|) tries to send A to at a rate proportional to 
ile + 7j)/2, but this effect may be compensated by the second term of the 
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RHS of fllSp . This means that nonzero solutions of the BO amphtude may 
be present below a critical value of (7e + 7^), i.e., when the scattering time 
is long enough. 

2.2. Boundary conditions 

Equations ffT7|) - fl2U]) must be solved together with the following boundary 
conditions at the cathode (x = 0): 

^^^ + aoF(0,t) = (J),(t), 

A{0,t) = 0, (32) 



and at the anode (x = L): 

dF{L,t) 



aLn{L,t)F{L,t) = {J),{t), 



dt 

A(L,t) = 0, (33) 



together with the voltage bias integral constraint: 

Y I F{x,t)dx = (p, (34) 
^ Jo 

where cxo and ctl are the dimensionless conductivities at the cathode and 
anode respectively, L is the nondimensional SL length and (pL is the nondi- 
mensional constant applied voltage bias. 

The contact conductivity at the cathode do must be selected so that aoF 
intersects the second branch of the drift velocity {J)g (F) depicted in Fig. [1] 
{J)q (F) is found by solving ( TT7|) for F, provided n = 1 and F is independent 
of X and t: 

i^Mn-j^. (35) 

This is a typical boundary condition that yields Gunn t ype self-sustained 
oscillations of the current in drift-diffusion SL models 17, 18, 22. [isj. 
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Figure 1: 0-averaged current density versus electric field for a spatially homogeneous 
stationary state at different temperatures. 



2.3. Initial conditions 

We select spatially uniform A and F as initial conditions: 

A(x,0) = ^^i^, (36) 



F(x,O) = 0, (37) 



where Ij{x) is the modified Bessel function The initial value of A is 

calculated by using fH5|) assuming that the initial values of u^^\ 13^^^ and 
6 are 0, A/(2A;bTo) and 0, respectively. The initial nondimensional mean 
energy is the lattice mean energy density: 

Eo = ^4^. (38) 



2kBTo 



Although the initial-boundary value problem (IBVP) to be solved is one 
dimensional, we need efficient numerical methods because it is very computa- 
tionally intensive: we need a large integration time to observe the Gunn-type 
self-oscillations of the current and the number of operations needed to com- 
pute the double Fourier coefficients /^o°\ f2-h fs-i ^^"^ f2-i given by 
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(125]) - fl26p is large. The numerical computation of the multipliers u and /3 
(which depend on both fast and slow time scales) is one of the bottlenecks of 
this problem. To calculate f2o'\ f2^i and f^^i in flTTl) and ffTS]) . we need 
to know u^^^ and /3'-°'' given by solving f l29|) . To calculate f2-i we need m*-^^ 
and given by dSU]). 



3. Numerical solution 



After replacing n = 1 + dF/dx from ( |20|) . the drift-diffusion equation ( !T7|) 
can be written in the following way: 
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The influence of the amplitude A in equation ([T 
the term f^o'^ in (H3!) contains it. On the other hand, in equation (fTS!) we 
can separate the diffusion term from the rest of the terms. Note that the 
amplitude A appears implicitly in the double Fourier coefficients , f^-i 

and f2-i- Therefore, we can write f lTS]) as follows: 
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For the numerical computation of u'^^^ and (3'^'^^ we can write equation (|29 
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in which 
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where 



K^= e^^°^''°'''cosh{u^^^k)dk (47) 
Jo 

Ks= j e^^"'' smh{u^^'^ k) sin kdk (48) 
Jo 

= = e'^*"' ^^'^ cosh(M(o)A;) cos k dk (49) 

For an efficient numerical computation of integrals P7|) -f l49|) we use an expan- 



sion as series of modified Bessel functions described in [Appendix A[ Once 



we have obtained u^^^ and (3^^^ from (H^ . we can get u^^^ and f3^^^ as shown 



in Appendix Appendix A 



3.1. Numerical scheme 

We have used an implicit numerical scheme to solve the partial differential 
equations f lT7|) and f lTSj) . In order to avoid numerical instabilities, our scheme 
employs a fixed point iteration method to find the numerical solution for F, 
A, (Jg) and J at each time step. 

3.1.1. Drift- diffusion equation for F 

To solve equation f|T7|) , we use a scheme similar to the one described 
and proved to converge in |2J] for a related problem involving partial differ- 
ential equations with an integral constraint. We use central differences for 
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approximating spatial derivatives, and the resulting differential equation is 
integrated in time by a first order implicit Euler method. This procedure 
leads to a system of + 2 linear equations for the + 1 values of the electric 
field F° f» F{jh, nr), j = 0, 1, . . . , A^, with the subscript j referring to space 
and the superscript n to time, plus {J)g. In order to save computational ef- 
fort we will set up the finite difference system of equations with a tridiagonal 
coefficients matrix, in the following way: 

a,F-_Y + + c.F^^l + d, = g„ j = l,...,N-l (50) 

The coefficients of fISU]) are: 

a, = -^Ar' + rBr' (51) 
hr 



g^ = h^Ff + h''TVf\ (52) 

where h = Ax = L/N, r = At and all the coefficients A, B, C and V of (139!) 
are evaluated at time 

The voltage bias integral constraint is solved by the composite Simpson's 
rule: 

Fo"+i + 4F^+' + 2F^+' + ... + 2F^+_\ + 4F;^-'i + F^+^ = 3<j)L/h. (53) 

The boundary condition at the injector contact is: 

{l + aoT)Fr'-r{jC=F^, (54) 

and at the collector contact: 

{h + aLT{h + - F;^+\))F;^+^ - rh ( J)^+^ = hF^^. (55) 

This system of A^ + 2 linear equations can be reduced to a simpler and smaller 
system, with the objective of finding a tridiagonal matrix, in the following 
way: 
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(J)g can be calculated directly from the boundary condition at the 
injector contact: 



(J) 
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(56) 



The field at the anode can also be expressed in terms of Fq^^ and F^_\: 
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We can make the following factorization of the system of linear equa- 
tions: 
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and vectors F, v, g and u are: 
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System can be efficiently solved by means of the following 

system with the same tridiagonal matrix T: 



Ty = g 
Tz = V 

After calculating y and z, we can obtain F^-^\ F and (J)^^^ 

u ■ y — 3(f)L/h + K,2 
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For each 9 G [— vr, vr], the nondimensional multipliers f3j~^^'^ and 

are obtained by solving f l46|) using the Newton-Raphson method with 

the Jacobian matrix f lA.ip . The Boltzmann distribution function Fourier 

coefficients ^/^o°^ j , (^f2,-i^ j are calculated from 

( 123|) -( 12^ using the composite Simpson's rule for all the integrals over 
k and over 6. 



3.1.2. Pamholic equation for A. 

It is important to find an accurate finite difference scheme to solve f llSp 
because the restitution parameters 7e,7j are close to their critical values; 
thus we will use a second-order implicit scheme. We approximate the time 
derivative of the complex amplitude A at time n + 1 as: 



dA 



n+l^ 



dA 



n+l 



dt 



n+l 



r 



^ /in 



(66) 



which is equivalent to a second-order, implicit, Runge-Kutta method (trape- 
zoidal rule). We use central difference for the space derivative in the second 
term of the RHS of (HH) clS cl whole block: 



5s 



d (fin 
dx\ l + iF 



(1) 



1 

2h 



n+l 



i+1 



B{0) 
2,-1 



6s 



(1) 

2,-1 



n+l- 



J-1 



1 + tF^:,' 



1 + 



(67) 
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Therefore, the numerical scheme for A will be as follows: 



T 



8ih 



/frB{0) . r (1) ffB{0) , r (1) 

[f2,-i + ^ 4,-1 ) .^^ [f2,-i + ^ rs'^ '_i ) \ 



8{1 + iF)h^ 32(l + iF)2/i2 ^ ) 

where j = l,...,A^ — 1. The boundary conditions at the contacts are: 

A^' = 0, (69) 
A'^^' = 0. (70) 

3.1.3. Algorithm 

1. For each time step t"^^^ do: 

(a) While the fixed point iteration does not converge: 
• For each node j = 1 . . . — 1: 

— For each 6 G [— 7r,7r]: calculate the Jacobian matrix el- 
ements from (lA.ip using the Bessel functions series de- 



scribed in Appendix Appendix A for the integrals P7|) - 
(glD and (jA3)-(Hl), and then obtain (/3W)°+^'^ and {u^^^Y-^^'' 
by the Newton- Raphson method from f j45l) . 

Calculate {fTY^'^ {f^-lY^' and (r^^li)"^' from dH, 

(1241) and (pT]) . The integrals over k and ^ are calculated 
by the composite Simpson's rule, with the previous values 
obtained of (/3(o));+i'^ and (^(0));+'''. 

Calculate coefficients A]+\ B]+\ and V]+^ from 

m-m. 



— Calculate coefficients Uj, bj, cj and gj from fl5Tl) - fl52l) . 

(b) Obtain -F""*"^ and (^7)^^^ from fl^ - flBS]) . the complex amplitude 
A""*"^ from fl68|) - fl70|) and the current density J from f|T9l) . 

(c) If the fixed point iteration does not converge, then go to step (a). 
2. Go to step (1). 
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3.2. Numerical results 

We have solved the hydro dynamic equations using parameter values sim- 



ilar to those for the superlattice with A = 16 meV in Ref. [25[. Then 
5 = 0.0053 and other parameter values are as in Table [H To obtain un- 
damped BOs, we have used •jej = 1.1269 so that (7e + 7j)/2 < 7crit- We 
consider a 50-period dc voltage biased GaAs-AlAs SL with lattice tempera- 
ture 70 K and dimensionless contact conductivities ctq = 1 and ai = 0.25. 
Initially, the mean energy density is Eq = 0.5501 and the profiles of A and 
F are uniform, taking on values of 0.5501 and 0.05, respectively. 



For a voltage bias = 0.05 {V = 0.166 V) and after a short transient 
that depends on the initial conditions, we observe coexisting BOs of frequency 
0.36 THz and Gunn type oscillations of frequency 11 GHz. BOs are stable 
because (7e + 7j)/2 < 7crit and Gunn type oscillations are a consequence of 
the periodic recycling and motion of electric field pulses from the cathode to 
the anode. Figure |2] shows several snapshots of the field and |y4| profiles of the 
Gunn type oscillation, and Fig. [3] exhibits the corresponding current density 
profiles for 6 = 0. While the amplitude of Gunn-type current oscillation is 
about 0.03 in nondimensional units (as seen in Fig. [2]^a) for the total current 
density averaged over the BOs), the BO part of the current oscillation has 
a larger amplitude of about 0.5 (as shown in Fig. |2]^c) for the modulus of 
A). Figure m illustrates the total current density ( IT^ of the coexisting 0.36 
THz Bloch and 11 GHz Gunn type oscillations, respectively. For each lattice 
temperature, there is a critical curve in the plane of restitution coefficients 
such that, for (7e + 7j)/2 > 7crit, BOs disappear after a relaxation time but 
they persist for smaller values of (7e + 7j) [11]. 

Figure O shows the profiles of F and 1^41 and Fig. [HI depicts the total 
current density at temperature 300 K, with Eq = 0.1529, for the same values 
of 7e,j and the other parameters. We find BOs but not the slower Gunn 
type oscillations. Whether Bloch and Gunn type oscillations coexist depends 
on the relative size of the diffusion (HT!) and convection (140!) terms in (139!) 
which, in turn, are controlled by the lattice temperature according to 
There is a critical temperature below which diffusion terms are sufficiently 
small compared to convective terms in f|T7j) and the electric field pulses are 
then periodically recycled when they reach the anode, originating the Gunn 
type oscillations. For larger temperatures, Bloch and Gunn type oscillations 
cannot occur simultaneously: When the electric field pulse reaches the anode, 
it remains stuck there and the electric field profile becomes the stationary 
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Figure 2: (a) 6'-averaged total current density vs time during coexisting Bloch and Gunn 
type oscillations at 70 K. (b) Field profile vs space at the times ti to <4 marked in (a), (c) 
Same for the complex BO amplitude profile. To transform the magnitudes in this figure 




Figure 3: Profiles of tlie nondimensional electron current density J„ and the modulus of 
the amplitude for BOs a.t 9 — corresponding to the instants ti — of Fig. [2][a). 



20 




Figure 4: (a) Total current density vs time during coexisting Bloch and Gunn type os- 
cillations at 70 K. (b) Fourier transform of the total current density showing two peaks 
corresponding to coexisting Bloch (0.36 THz) and Gunn type (11 GHz) oscillations. 



state shown in Fig. |5l Note that the largest peak in the current spectrum in 
Fig. El^b) occurs at a lower frequency (0.26 THz) than in the case of lower 
lattice temperature shown in Fig. Hl^b). 

3.3. Spurious convection 

Had we used first-order backward differences for the space derivative in 
( IT8l) instead of second-order ones (1671) . the computing time to get comparable 
results would have had to increase significantly: unless we use four times 
smaller values of Ax, the scheme with first-order differences gives rise to 
spurious convection terms in A. Fig. [7] shows the spatial profile of \A\ at a 
given time during BOs, calculated using first-order backward differences in 
( 1T8|) for two values of (7e + 7j) < (7e + 1j)crit- The dashed lines in Figs. W^a) 
and (b) show that BOs extend to the whole SL for sufficiently small values 
of (7e + 7j), no matter the step size. However, when (7e + 7j) approaches 
the critical value from below (solid line), the BOs are confined to part of 
the SL as in Fig. W^a) unless the step size is sufficiently small as in Fig. 
l^b). Fig. W^a) shows that, for larger step size, there appears a spurious 
convection in A that extends the region where A = from the cathode and 
it confines the BOs to the region closer to the receiving contact (anode). The 
spatial interval where A = moves from cathode towards the anode as time 
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(a) (b) 




0.5 1 1.5 2 0.5 1 1.5 2 

x/[x] x/[x] 



Figure 5: (a) Electric field profile vs space for the stationary state at 300K. (b) Same for 
the modulus of the BO complex amplitude. To transform the magnitudes in this figure to 
dimensional units, use Table [1] 




Figure 6: (a) Total current density vs time during Bloch oscillations at 300K. (b) Fourier 
transform of the total current density showing only one peak corresponding to BOs (0.26 
THz). The zero-frequency constant corresponding to the time average of the total current 
density has been subtracted. 
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Figure 7: Modulus of the nondimensional complex BO amplitude at a given time, calcu- 
lated using first-order backward differences to solve ([T^ with different space steps. 



increases. The phenomenon of the spurious convection occurs for restitution 
coefficients close to their critical values as in Ref. |10j]. 

4. Convergence of the method 

We have verified the convergence of the method in terms of the spatial and 
temporal mesh size Ax and At, by means of numerical simulations. Figure [H] 
shows the evolution of {J)g for different values of Ax and At. Convergence 
to the solution is more sensitive to the size of the space step Ax than to the 
size of the time step At. We observe in Fig. IHt^a) that if we use two or more 
integration nodes per SL period (i.e., Ax < 0.022), the difference between 
the numerical solutions is small and the artifact that appears for large space 
step. Ax = 0.044 is eliminated for smaller step sizes. On the other hand. Fig. 
IHl^b) shows that decreasing the time step by half has little effect in improving 
the convergence of the scheme: the artifact for large Ax = 0.044 still appears 
whereas for smaller space steps, the improvement provided by a smaller time 
step is quite small. 

We also analyzed the average number of fixed point iterations necessary 
for convergence of the numerical scheme at each time step and checked that 
the fixed point iteration is contractive. In the case where the solution consists 
of Gunn-like self-oscillations, only two or three iterations are needed for each 
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Figure 8: 0-averaged current density versus time for different values of space and time 
steps. In (a) we use At = 0.01. 
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Figure 9: Error of 0-averaged current density (in the L^-norni) for different values of Aa; 
and At. In (a) we use At 0.01. In (b) we use Ax = 0.011. 
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time step during most of one period. However when the field pulses reach 
the SL end and new pulses are created at the cathode, more than three 
iterations are needed to attain a reliable result. We have observed that the 
number of iterations depends directly on the size of the time step At. Table 
[2] summarizes the main results of the convergence analysis for a simulation 
of t = 200 [t] = 376 ps. It can be seen that for larger At, more fixed point 
iterations are needed at each time step, therefore the total computation time, 
for a given Ax, does not depend too much on At as long as it is small enough 
for the convergence of the iterative scheme. 

Figure [9] demonstrates the convergence of our numerical scheme. Figures 
M (a) and (b) depict the time- averaged error of {J)g (in the L^-norm) as a 
function of the spatial and time discretizations, respectively. We observe 
that the convergence order is approximately quadratic 1.88) in space and 

1.92) in time. 



N 


Ax 


At 


# iterations/^ steps (t) 


Total 7^ iterations 


Mfiops/^^ iterations 


50 


0.044 


0.03 


3.011 


20073 


0.554 


50 


0.044 


0.02 


2.123 


21230 


0.554 


100 


0.022 


0.03 


2.987 


19913 


0.570 N 


100 


0.022 


0.02 


2.138 


21380 


0.570 N 


150 


0.015 


0.03 


2.945 


19633 


0.572 


150 


0.015 


0.02 


2.153 


21530 


0.572 


200 


0.011 


0.03 


2.913 


19420 


0.575 N 


200 


0.011 


0.02 


2.334 


23340 


0.575 A^ 



Table 2: Average number of fixed point iterations and Mflops. 



We can also observe in Table |2] that the average number of fioating point 
operations is almost proportional to the number of space integration nodes 
(A^). Considering the average number of iterations per time step, the total 
computational cost for a simulation with A^ = 200, At = 0.02 and Ngteps = 
10000 (Cost ^ 0.575 x 10^ NNueraUonsN steps) is of the order of 2.4 x 10^^ 
fiops. The computation time in a computer with Intel(R) Core(TM)2 Duo 
CPU E6750 @ 2.66 GHz processor is about 30 hours. These estimates show 
that although the problem is one dimensional, it is rather computationally 
intensive and therefore it is important to optimize the numerical algorithm. 
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5. Conclusions 

We have presented a finite difference method to numerically solve a new- 
type of hydrodynamic equations that arise in the theory of spatially inhomo- 
geneous Bloch oscillations in semiconductor SLs. The hydrodynamic equa- 
tions describe the evolution of the electric field, electron density and the 
complex envelope of the Bloch oscillations for the electron current density 
and the mean energy density. They contain averages over the Bloch phase 
which are integrals of the unknown electric field. These equations are de- 
rived by singular perturbation methods from a Boltzmann-Poisson transport 
model of miniband SLs with inelastic collisions. Among the solutions of the 
hydrodynamic equations at low lattice temperature, there are spatially inho- 
mogeneous Bloch oscillations coexisting with moving electric field domains 
and Gunn-type oscillations of the current. The latter oscillations disappear 
at higher temperatures (300 K). These novel Bloch-oscillation solutions are 
found for restitution coefficients in a narrow interval below their critical val- 
ues and disappear for larger values of the restitution coefficients. 

To solve the averaged hydrodynamic equation for the time evolution of the 
complex BO amplitude, we use an efficient implicit second-order numerical 
scheme that uses a fixed-point iteration process to avoid numerical instabili- 
ties. In the case of the drift-diffusion equation for the electric field, only one 
tridiagonal matrix needs to be inverted to solve the implicit scheme, which 
results in a greatly decreased computational time. Double integrals entering 
the averaged hydrodynamic equations are calculated by means of expansions 
in modified Bessel functions. We use numerical simulations to ascertain the 
convergence of the method. If the complex amplitude equation is solved us- 
ing a first order scheme for restitution coefficients near their critical values, 
a spurious convection arises that annihilates the complex amplitude in the 
part of the superlattice that is closer to the cathode. This numerical artifact 
disappears if the space step is appropriately reduced or we use the second- 
order numerical scheme. 
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suggesting the use of expansions in series of Bessel functions. 
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Appendix A. Numerical method for obtaining u and /3 



Solving f l45p by means of the Newton- Raphson method with the following 
Jacobian matrix: 



J 



a ( £c£ 

d ( Ks 



9/3(0) I K, 



a 



Ks 

9/3(0) 



(A.l) 



we obtain the numerical values of u*^"-* and Thus we need to calculate 
numerically the integrals fH7j) - fH9|) . as well as the following ones: 

^'^^ 3^'°' '-^'^ cosHu^^h) cos^ k dk, 



K. 



CI3I3 



K. 



g/3(o)cosfc^gi^j^^^(o)^^ cos/trfA;, 







,/3(0) cos A: 



OK, 



^/3(0) cosfc 



3/3(0) c;osA: 



sinh('u'-°-*fc) sin k cos k dk, 
k cosh.{u^'^^ k) sin k dk, 
k sinh(M(°''fc) dk. 



(A.2) 
(A.3) 
(A.4) 
(A.5) 
(A.6) 



We also need the integrals dSD-dinD and ([0|) - ([0|) to obtain m^^) and (3^ 
from flA.lOl) and flA.lip . and therefore an efficient method is required for 
their numerical computation. If we expand e^^°'^°'^'^ as a series of modified 
Bessel functions, these integrals can be expressed as the following series (we 
omit the argument (3^^^ of the Bessel functions): 



K, 



sinh(7ru(0-') 

K, = sinh(7ru(°)) 
sinh(7rM(o) 



/o + 2(«(°))^5^(-l) 



I. 



1 



(o)^2 



f + (u(0))2 

1+J 



K. 



c/3 



K, 



sinh(7r-u(0'' 



i=i 

00 



1-J 



(j + 1)2 + (m(0))2 (1-j)2 + (^(0))2 



^ f + (m(0))2 



(0) 



2(.(o))^$: 



^ f + (m(0))2 
i=l J V / 
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K. 



cup 



i=i 



y(0) Z^V (^-2 + (^(0))2)2 \ J 



2:0 

/3(o) 



Ksp = smh(7rM(°)) 



^ j2 + {u(^)f 



1 + 



1+J 



+ 



1-J 



su — ^0 



7r(l + cosh(7ruW) - 2mW sinh(7ru(o)) 
In : — tt^ttittt: 



(1+ (m(0))2)^ 



1+J 



(j + 1)2 + {u(^)y 



1-J 



(1-j)' + (m(°))V " 

= po))^ 1 ^0 + 2(m^ ^) 2^(-ir-^^ 



j2 + (m(0))2 



+ ^^E(-l) 



(0) 



(j2+ (m(0))2)2' 



where Ij = Ij{(3^^^) is the modified Bessel function of the first kind with 
index j and argument (3^^\ Since \u^^^\ and \(3^^^\ are less than 3.5, with only 
13 terms of the previous Bessel series we can have an error less than 10~^ 
percent. 

The numerical value of m*^^-* and can be obtained from (130|) considering 
that /^(i) is 

rW=7™(Gi/3« + G2««), (A.7) 

in which Gi and G2 are: 



(A.8) 
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\Kc ) 

From flA.7p we can derive j^^^^^ ^'^'^ then u^^^ and /^'-^^ can be explicitly 
obtained by solving a linear system of equations: 

(1) ^ K {QuR2 - Q21R1) . 
n (Q11Q22-Q12Q21)' ^ ■ ' 

n(i) _ (Q22 ^1 - Q12 R2) (^^^\ 
' n {QuQ22-Qi2Q2iy ^ ■ ^ 

where Ri and R2 are the real and imaginary parts of the RHS of fl30|) and 
the Qij are: 

Qll = KcKclijB — {Kcpf , 
Q12 = Kc Kcul3 ~ Kci3, 

Q21 = —Kc Ksi3 + Ks Kcp, 
Q22 = —Kc Ksu + Kg Kcu- 
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